Sudan_catch_analysis

Analysis of Sudan survey data 2012 - 2013

This R-markdown document contains all r-code used to carry out the analysis for the paper on fish catch rates and species densities (species richness) along the Red Sea coast of Sudan. The analysis is based on three surveys: November 2012, May 2013 and November 2013.

Here the analysis is structured around the 7 management regions of the Sudanese Red Sea coast as follows:

• Present the region: Bathymetric map
• Sampling locations
• Table of effort (diffent gear types per management areas.)
• Depth distribution of traps by survey and area
• Catch CPUE by area and survey
• Catch CPUE distribution by gear, by species groups (should perhaps do the same for trophic groups)
• CPUE by trophic group pr area and survey
• CPUE by depth and tropich group and survey
• Biodiversity /species densities
• Statistical analysis
• test of normality shows that we need to use non-parametric and zero-adjusted approaches.
• CPUE as a function of area, depth, survey, gear (Zero-adjusted model)

All code and data are stored on GitHub

Libraries, data etc.

Libraries, dependencies and functions

Load and manipulate data

Reading catch, station and traits data.

Neither station data, nor catch data has a complete depth record, but by combining depth para from each we can get a complete depth parameter for all stations.

Sudan map and management areas

Loads the map data for Sudan, loads the management areas from shape-file etc.

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/eriko/Documents/GitHub/Sudan2019/Sudan-master/sudan_management_areas", layer: "sudan_regions"
## with 7 features
## It has 3 fields
## Integer64 fields read as strings:  id

Allocating catch positions to management areas

Adds ‘Area’ to each line in the catch table.

Modifying dataset and creating subset for stations with catches

Adding traits to ‘catch’ data set (without 0-catch stations)

Some more data wrangling that adds the traits from the traits table to the catch data, using only the station with catches (there are no traits for ‘NOCATCH’ species).

Also, only select traps, gillnets and handlines as these were the only gear with sufficient numbers and consistent use to be analyzed.

Lastly adds number of gear deployed at each station to each line.

Bathymetric map of Sudan

Fig. 1 in MS

A nice, bathymetric map of Sudan with management areas overlaid.

Currently fonts and colours are adjusted for use on the ICES ASC 2019 poster, removing the bathy-contour lines (se ‘lwd’ to 0), and varying the colours of management area text to white for areas 1-4.

## quartz_off_screen 
##                 2

Stations plotted on maps, pr survey

Faceted map plotting position of all catch stations for each of the three surveys.

MApoints<-data.frame(lon=c(37.2, 37.53, 37.5, 37.4, 37.3, 37.9, 38.2), lat=c(21.7, 20.92, 20.3, 19.9, 19, 19.3, 18.5), area=c(1,2,3,4,5,6,7))

sudan.map<-ggplot(catch, aes(x=lon, y=lat, group=survey_m)) + geom_point(shape="+", size=6, aes(group=survey_m, colour= survey_m)) + facet_wrap(~survey_m, ncol=3) + geom_path(data=ManageAreas.f,colour="dodgerblue3", aes(x=long, y=lat, group=id)) + geom_polygon(data=world, colour="gray35", fill="gray85", aes(x=long, y=lat, group=survey_m)) +  coord_cartesian(xlim = c(36.5, 39), ylim=c(17.7, 22.5)) + theme_classic() + geom_point(data=ps, size=2, colour="gray35", aes(x=x, y=y)) + geom_text(data=ps, label="PZU", size=3,  hjust=1, vjust=-1.2,  aes(x=x, y=y, group=survey_m)) +geom_text(data=MApoints, aes(x=lon, y=lat, group=area, label=area)) +scale_colour_brewer(type = "seq", palette = "Dark2", name="Survey month & year", labels=c("Nov. 2012", "May. 2013", "Nov. 2013")) + theme(legend.position="none")
sudan.map

ggsave("./figs/station_map.tiff")

## Saving 7 x 5 in image

Station information

Station table

Table describing the sampling effort in each area pr survey, number of different gear types, max, min and average depth, number of traps with / without catches.

#' TABLE OF SAMPLING EFFORT PR. MANAGEMENT AREA
#' (new for resubmission to PlosOne)
#' --------------------------------------------

s <- unique(catch.unmod$survey)
a <- as.factor(1:7)
ay_info <- data.frame(survey=integer(), id=integer(), Ntraps=integer(), Nhl=integer(), NGn=integer(), TBhrs=double(), HLhrs=double(), GNhrs=double(), DepthAvg=double(), DepthMax=double(), DepthMin=double()) 
cn <- colnames(ay_info)
# depth for TB only (GN are at the surface, and HL not measured)

for (i in 1:length(s)) { 
  y <- subset(catch.unmod, survey==s[i])
  for (j in 1: length(a)){
    ay_1 <- data.frame(survey=integer(), id=integer(), Ntraps=integer(), Nhl=integer(), NGn=integer(), TBhrs=double(), HLhrs=double(), GNhrs=double(), DepthAvg=double(), DepthMax=double(), DepthMin=double())
    ya <- subset(y, id==a[j])
    ay_1 <- rbind(ay_1, c(
      as.integer(as.character(ya$survey[1])), 
      as.integer(as.character(ya$id[1])), 
      length(unique(subset(ya, gear=="TB")$station)),  
      length(unique(subset(ya, gear=="HL")$station)), 
      length(unique(subset(ya, gear=="GN")$station)), 
      sum(subset(ya, gear=="TB")[!duplicated(ya$station),]$Fhrs, na.rm=TRUE),  
      sum(subset(ya, gear=="HL")[!duplicated(ya$station),]$Fhrs, na.rm=TRUE), 
      sum(subset(ya, gear=="GN")[!duplicated(ya$station),]$Fhrs, na.rm=TRUE), 
      mean(subset(ya, gear=="TB")[!duplicated(ya$station),]$depth, na.rm=TRUE), 
      max(subset(ya, gear=="TB")$depth, na.rm=TRUE), 
      min(subset(ya, gear=="TB")$depth, na.rm=TRUE)    ))
    colnames(ay_1) <- cn
    ay_info <- rbind(ay_info, ay_1)
    colnames(ay_info) <- cn
  }
}

formattable(ay_info)

survey	id	Ntraps	Nhl	NGn	TBhrs	HLhrs	GNhrs	DepthAvg	DepthMax	DepthMin
2012901	1	22	0	0	694.066	0.000	0.000	42.45455	142	13
2012901	2	54	3	3	721.549	78.000	62.000	40.89744	71	0
2012901	3	26	0	1	451.183	0.000	14.000	30.80769	95	8
2012901	4	5	0	0	77.084	0.000	0.000	22.60000	54	10
2012901	5	31	0	4	677.896	0.000	78.317	21.32258	30	7
2012901	6	36	0	1	850.417	0.000	38.250	32.38889	88	0
2012901	7	31	0	8	712.200	0.000	135.849	31.06667	66	5
2013002	1	29	0	1	420.163	0.000	24.000	31.48148	70	5
2013002	2	81	3	5	1102.379	377.500	221.569	27.09091	145	5
2013002	3	32	0	1	545.811	0.000	24.000	28.87500	60	0
2013002	4	13	0	10	160.266	0.000	137.815	29.50000	67	9
2013002	5	33	2	5	208.899	192.000	195.766	20.46154	50	7
2013002	6	45	1	2	642.413	156.000	78.000	34.70270	88	9
2013002	7	39	2	0	666.410	186.000	0.000	33.46154	76	5
2013005	1	23	1	4	271.551	3.000	84.000	29.82353	80	10
2013005	2	57	2	6	500.101	111.000	156.000	38.27273	80	7
2013005	3	9	0	2	123.099	0.000	36.000	26.50000	70	9
2013005	4	16	2	4	151.184	4.500	142.767	32.90000	68	12
2013005	5	30	2	3	317.498	7.000	146.947	25.52381	65	11
2013005	6	40	4	2	443.983	31.501	71.915	40.14815	89	6
2013005	7	22	0	2	171.784	0.000	203.171	34.75000	54	11

write.csv2(ay_info, "./figs/Area_Station_table.csv")

Number of set traps per depth

Violin plot of depths at wich traps were set. I prefer this visualization of depth ranges of the traps, rather than a (boring) table.

c2 <- catch %>% group_by(ss) %>% summarise(first(depth), first(gear), first(id), first(survey), first(survey_m))
colnames(c2) <- c("ss", "depth", "gear", "id", "survey", "survey_m")

dp <- ggplot(subset(c2, gear=="TB" & depth!=0), aes(as.factor(id), depth))

dp <- dp +  geom_violin(scale="width", aes(fill=as.factor(id)), draw_quantiles = c(0.5)) + scale_y_reverse() + labs(x="Management areas", y="Depth (m)") + scale_fill_manual(name="Management area",values=sequential_hcl(n = 7, h = c(0, -100), c = c(80, NA, 40), l = c(53, 75), power = c(1, 1), rev = TRUE, register = ), guide=FALSE ) +theme_bw()  +facet_wrap(~as.factor(survey_m))

dp

tiff(file="./figs/trap_area_depth.tiff", width=3000, height=1500, res=400, pointsize=10, compression=c("none"))
dp
dev.off()

## quartz_off_screen 
##                 2

Analysis of catches by traits and other factors

Catch rates pr area

#´ Explore: Trap Catch rates (CPUE) by depth and survey
#' not used in MS
#d3 <- ggplot(subset(catch, gear=="TB"), aes(depth, CPUEw))
#d3 + geom_point() + geom_smooth(method="auto")  +facet_wrap(~as.factor(survey_m))

#' main catch by depth plot pr area
ca.plot<-ggplot(subset(catch, gear=="TB" & depth!=0), aes(as.factor(id), CPUEw))

ca.plot <- ca.plot + geom_violin(scale="width", aes(fill=as.factor(id)), draw_quantiles = c(0.5))  + labs(x="Management areas", y="CPUE kg/hrs fishing") + scale_fill_manual(name="Management area",values=sequential_hcl(n = 7, h = c(0, -100), c = c(80, NA, 40), l = c(53, 75), power = c(1, 1), rev = TRUE, register = ), guide=FALSE ) +theme_bw() +facet_wrap(~as.factor(survey_m))
ca.plot

tiff(file="./figs/catch_rates_traps_area_survey.tiff", width=3000, height=2000, res=400, pointsize=10, compression=c("none"))
ca.plot
dev.off()

## quartz_off_screen 
##                 2

CPUE by gear type - by family and trophic group

# by genus (families)
catch.3gear<-subset(catch.unmod, gear=="GN" | gear=="HL" | gear =="TB")

gear.plot <- ggplot(catch.3gear, aes(capwords(tolower(as.character(catch.3gear$FamGroup))), fill=survey_m)) + geom_bar(binwidt=0.2, aes( weight=weight)) +theme_bw() + facet_grid(vars(gear), vars(survey_m)) + coord_flip() +xlab("") + ylab("kg")  + guides(fill=guide_legend(title="Survey")) + theme(legend.position = "none")

## Warning: Ignoring unknown parameters: binwidt

gear.plot

ggsave("./figs/gearplot_species.tiff")

## Saving 7 x 5 in image

# by trophic group
gear.plot2 <- ggplot(subset(catch.traits, TGShort!="NA"), aes(TGShort, fill=survey_m)) + geom_bar(binwidt=0.2, aes( weight=weight)) +theme_bw() + facet_grid(vars(gear), vars(survey_m)) + coord_flip() +xlab("") + ylab("kg")  + guides(fill=guide_legend(title="Survey")) + theme(legend.position = "none")

## Warning: Ignoring unknown parameters: binwidt

gear.plot2

ggsave("./figs/gearplot_trophic_group.tiff")

## Saving 7 x 5 in image

CPUE by trophic group by area & by survey (year)

FIG 5 in Feb 2019 version of MS

Presents the CPUE pr hr, for each management area and survey, standardized by the number of traps set in each area X survey combination.

Trap catches were dominated by carnivores, occuring in all area X survey combinations, followed by invertivores who also occurred everywhere, albeight in lower densities than the carnivores.

Area 2, 5 and 6 seem to have the highest diversity of trophic groups, but this has to be calculated as functional diversity, using the FD package (or similar).

trophic.plot <- ggplot(subset(catch.traits, TGShort!="NA" & gear=="TB"), aes(TGShort, CPUEw/NArea))

trophic.plot <- trophic.plot + geom_col(aes(fill=TGShort)) + facet_grid(vars(id), vars(survey_m)) + scale_fill_manual(name="Trophic group",values=qualitative_hcl(n = 5, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x=" ", y="CPUE (kg/hrs fishing) / (number of traps pr area)") + theme_bw() 
trophic.plot

tiff(file="./figs/trophic_group_cpue.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
trophic.plot
dev.off()

## quartz_off_screen 
##                 2

Trophic group CPUE by depth

A modified version of Figure 6 in Feb 2019 MS.

Basically a more complex version of the previous plot. Adds the depth dimension to the plot, and the surveys as colours.

Shows uniform catch rates of carnivores and invertivores from 0-50m.

All surveys seem to overlap, so no difference in trophich group depth dependent catch rates between the surveys.

depth.plot <- ggplot(subset(catch.traits, gear=="TB" & TGShort!="NA"), aes(depth, CPUEw)) 

depth.plot <- depth.plot + geom_point(aes(colour=survey_m))  + facet_grid(vars(id),vars(TGShort)) +  scale_colour_manual(name="Survey",values=qualitative_hcl(n = 3, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x="Depth (m)", y="CPUE (kg(hrs)")  +theme_bw()  

depth.plot

tiff(file="./figs/trophic_group_cpue_depth.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
depth.plot
dev.off()

## quartz_off_screen 
##                 2

Trophic level in Trap catches

TMean <- subset(catch.traits, gear=="TB" & TrophicLevel!="NA") %>% group_by(id, survey_m) %>% summarise(TM=mean(na.omit(TrophicLevel)))
CMean <- subset(catch.traits, gear=="TB" & TrophicLevel!="NA") %>% group_by(id, survey_m) %>% summarise(CM=mean(CPUEw))
TMean$CM <- CMean$CM

t0.plot <- ggplot(subset(catch.traits, gear=="TB" & TrophicLevel!="NA"), aes(TrophicLevel, CPUEw)) 

t0.plot <- t0.plot + geom_point(aes(colour=TGShort)) + geom_vline(data = TMean, aes(xintercept = TM)) + facet_grid(vars(id),vars(survey_m)) +  scale_colour_manual(name="Trophic group",values=qualitative_hcl(n = 5, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x="Trophic level", y="CPUE (kg(hrs)")  +theme_bw()  
t0.plot

tiff(file="./figs/trophic_level.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
t0.plot
dev.off()

## quartz_off_screen 
##                 2

Maxlength in catches

t0.plot <- ggplot(subset(catch.traits, gear=="TB" & TrophicLevel!="NA"), aes(MaxLength, CPUEw)) 

t0.plot <- t0.plot + geom_point(aes(colour=TGShort)) + facet_grid(vars(id),vars(survey_m)) +  scale_colour_manual(name="Trophic group",values=qualitative_hcl(n = 5, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x="Max lenght(cm)", y="CPUE (kg(hrs)")  +theme_bw()  
t0.plot

tiff(file="./figs/max_length.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
t0.plot
dev.off()

## quartz_off_screen 
##                 2

Place in water column

catch.traits$WaterCol[catch.traits$WaterCol=="Demersal"] <- "demersal"
tl.plot <- ggplot(subset(catch.traits, WaterCol!="NA" & gear=="TB"), aes(WaterCol, CPUEw/NArea))

tl.plot <- tl.plot + geom_col(aes(fill=WaterCol)) + facet_grid(vars(id), vars(survey_m)) + scale_fill_manual(name="Place in water column",values=qualitative_hcl(n = 6, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x=" ", y="CPUE (kg/hrs fishing) / number of traps pr area") + theme_bw() + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
tl.plot

tiff(file="./figs/water_col.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
tl.plot
dev.off()

## quartz_off_screen 
##                 2

Diel activity

t2.plot <- ggplot(subset(catch.traits, DielActivity!="#N/A" & gear=="TB"), aes(DielActivity, CPUEw/NArea))

t2.plot <- t2.plot + geom_col(aes(fill=DielActivity)) + facet_grid(vars(id), vars(survey_m)) + scale_fill_manual(name="Diel activity",values=qualitative_hcl(n = 6, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x=" ", y="CPUE (kg/hrs fishing) / number of traps pr area") + theme_bw() + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
t2.plot

tiff(file="./figs/diel_act.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
t2.plot
dev.off()

## quartz_off_screen 
##                 2

Gregariousness

t3.plot <- ggplot(subset(catch.traits, Gregariousness!="#N/A" & gear=="TB"), aes(Gregariousness, CPUEw/NArea))

t3.plot <- t3.plot + geom_col(aes(fill=Gregariousness)) + facet_grid(vars(id), vars(survey_m)) + scale_fill_manual(name="Gregariousness",values=qualitative_hcl(n = 6, h = c(0, 295), c = 80, l = 60, register = ) ) + labs(x=" ", y="CPUE (kg/hrs fishing) / number of traps pr area") + theme_bw() + 
  theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
t3.plot

tiff(file="./figs/gregariousness.tiff", width=3000, height=3000, res=400, pointsize=10, compression=c("none"))
t3.plot
dev.off()

## quartz_off_screen 
##                 2

---

Biodiversity (species densities, trophic diversity etc)

Aggregated Species density per area for (traps only)

Species per area pr hour of fishing, traps only.

The highest number of species was observed in area 2 (36), followed by area 3, 6 and 7 (25 in each).

Species density (number of species pr hour of fishing - soak time of traps) ranged from 0.008 (areas 1, 2 and 6) to 0.019 (area 4). Area 4 had the highest species density across all surveys (figure below), while the area with the lowest species density varied by survey: area 1 in November 2012 , area 7 in May 2013 , and area 2 in November 2013 .

sp.pr.area <- subset(catch.unmod, gear=="TB") %>% group_by(id, survey_m) %>% summarise(n_distinct(species), n_distinct(ss), sum(Fhrs))
colnames(sp.pr.area) <- c("id", "survey","Nsp","NTraps", "Fhrs")

sp.pr.area$Species_by_NoStations <- (sp.pr.area$Nsp-1) / sp.pr.area$NTraps
sp.pr.area$Species_Hrs <- (sp.pr.area$Nsp-1) / sp.pr.area$Fhrs
write.csv2(sp.pr.area, "./figs/species_density_area.csv")

#calculate species density pr area (all years combined)
sp.pr.area2 <- subset(catch.unmod, gear=="TB") %>% group_by(id) %>% summarise(n_distinct(species), n_distinct(ss), sum(Fhrs))
colnames(sp.pr.area2) <- c("id","Nsp","NTraps", "Fhrs")
sp.pr.area2$Species_by_NoStations <- (sp.pr.area2$Nsp-1) / sp.pr.area2$NTraps
sp.pr.area2$Species_Hrs <- (sp.pr.area2$Nsp-1) / sp.pr.area2$Fhrs
formattable(round(sp.pr.area2, digits=3))

id	Nsp	NTraps	Fhrs	Species_by_NoStations	Species_Hrs
1	18	74	2036.710	0.230	0.008
2	36	192	4483.859	0.182	0.008
3	25	67	1471.675	0.358	0.016
4	12	34	593.399	0.324	0.019
5	19	94	2097.324	0.191	0.009
6	25	121	3114.607	0.198	0.008
7	25	92	2485.890	0.261	0.010

#' Species Density plot
species.plot<-ggplot(sp.pr.area, aes(x=id, y=Species_Hrs)) + geom_bar(stat="identity", fill="steelblue") +theme_bw() + scale_x_reverse(breaks=c(1,2,3,4,5,6,7)) 
sphrs_unmod<- species.plot  + coord_flip() + theme(legend.title = element_text(size=14, face="bold")) + theme(legend.text = element_text(size=12)) + ylab("Number of species caught in traps / median hours of  fishing") + xlab("Management areas from N (1) to S (7)") + facet_wrap(~survey)
sphrs_unmod

ggsave("./figs/species_density.png")

## Saving 7 x 5 in image

Species density pr station - Zero Adjusted GAM model

The variability and relative difference in rank of species diversity per area as observed in the FIGURE ABOVE was confirmed by the zero-adjusted GAM model (AIC: -4161, df: 12) were all surveys were significant factors in the model. In addition, area 1, 2, 3 and 4 were also significant factors affecting species density, while the non-significant factors were: depth, and areas 5, 6 and 7.

spue.df <-  subset(catch.unmod, gear=="TB") %>% group_by(ss) %>% summarise(first(survey), first(id), n_distinct(species), first(depth), sum(Fhrs))
colnames(spue.df) <- c("id", "survey", "area", "Nsp", "depth", "Fhrs")
spue.df$Species_Hrs <- spue.df$Nsp / spue.df$Fhrs
#hist(spue.df$spue)

# zero-adjusted model:
sd.mod<-gamlss(Species_Hrs~depth+factor(area)+factor(survey),family=ZAGA, data=spue.df) # AIC: -4161

## GAMLSS-RS iteration 1: Global Deviance = -4185.753 
## GAMLSS-RS iteration 2: Global Deviance = -4185.753

# gamlss(Species_Hrs~factor(area)+factor(survey),family=ZAGA, data=spue.df)  # AIC: - 4163 - slightly better than the one above, but want to show the one with 'depth' as a factor. 
summary(sd.mod)

## ******************************************************************
## Family:  c("ZAGA", "Zero Adjusted GA") 
## 
## Call:  gamlss(formula = Species_Hrs ~ depth + factor(area) +  
##     factor(survey), family = ZAGA, data = spue.df) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                         Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)           -3.0456597  0.0282803 -107.695  < 2e-16 ***
## depth                  0.0003143  0.0003925    0.801 0.423536    
## factor(area)2          0.0952550  0.0257799    3.695 0.000238 ***
## factor(area)3          0.1083585  0.0319608    3.390 0.000740 ***
## factor(area)4          0.0834147  0.0390748    2.135 0.033148 *  
## factor(area)5          0.0540349  0.0296866    1.820 0.069184 .  
## factor(area)6         -0.0337542  0.0278053   -1.214 0.225200    
## factor(area)7         -0.0464743  0.0294473   -1.578 0.114993    
## factor(survey)2013002  0.2071281  0.0174699   11.856  < 2e-16 ***
## factor(survey)2013005  0.2470968  0.0190004   13.005  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.66955    0.02708  -61.66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    -16.2      129.2  -0.125      0.9
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  674 
## Degrees of Freedom for the fit:  12
##       Residual Deg. of Freedom:  662 
##                       at cycle:  2 
##  
## Global Deviance:     -4185.753 
##             AIC:     -4161.753 
##             SBC:     -4107.594 
## ******************************************************************

Functional diversity

The ‘dbFD’ function only works when run in-line, not the whole code-chunk.

We used the same approach as Stuart-Smith et al to exclude very few observations of a single species. and excluded all species observations occuring in 2 or fewer stations (traps).

To avoid dbFD crashing ‘m’ was set in the ‘dbFD’ call (the number of PoCA axes kept as traits to do the FRic analysis). Several levels of ‘m’ were tested, and this could probably be tested further (to higher levels of ‘m’)

Results for n_spec > 2, m=10

Results based on the output from the dbFD analysis:

No area came out with the highest or lowest FD across all four metrics, but some patterns emerged:

• Area 1 was highest in FEve and FDiv, but lowest in FDis and RaoQ
• Area 2 was consistently high in all metrics, and highest in RaoQ
• Area 3 was high inn all metrics, but never the highest
• Area 4 was lowest in RaoQ, low in FDis, medium in FDiv and FEve
• Area 5 was high in all metrics, but never the highest
• Area 6 was high in FEve, but the lowest in all other three
• Area 7 was lowest in FEve, low in FDis, medium in FDiv and RaoQ

RaoQ is the most used metric (the one used by Stuart-Smith et al), and this is the metric we will use and discuss in the MS. It supports the analysis of species densities for Area 2, although Area 6 has the lowest RaoQ which indicates that although the number of species is higher, their functional diversity is lower (or what we sampled was lower) compared to the other areas.

Stuart-Smith tests various RaoQs by taking out one trait at a time and this way ranking the traits according to R2.

#ss_list <- levels(catch1$survey)
sp_occur <- catch1  %>% group_by(Sci_name) %>% summarise(n())
colnames(sp_occur) <- c("Sci_name", "n_spec")

# wrangle traits into a traitstable that can be used in FD
missing <- setdiff( x=traits$x1, y=subset(sp_occur, n_spec>2)$Sci_name)
t5 <- traits[!traits$x1 %in% missing,]
tr <- t5[,c(4:10)]
rownames(tr) <- t5$x1
tr[(tr=="#N/A")] <- c(NA)
tr <- droplevels(tr)
tr$diel<-as.integer(tr$Diel.Activity)
tr$diel[(tr$diel==2)] <- c(0)
tr <- tr[,c(1,2,3,4,6,7,8)]

# wrangle 'catch1' to aggregate CPUE pr area (row) pr species (columns)
missing <- setdiff( x=catch1$Sci_name, y=rownames(tr))
c5 <- catch1[!catch1$Sci_name %in% missing,]
c_melt <- melt(c5[,c(28,11,22)], id=c("id", "Sci_name"))

#convert abundance to integer for better calculation performance
med_int1000 <- function(x) as.integer(median(x)*1000)
abund <- dcast(c_melt, id ~ Sci_name, value.var="value", fun=med_int1000) 
# to calculate with raw abundance data, use line below:
#abund <- dcast(c_melt, id ~ Sci_name, value.var="value", fun=median)
abund[is.na(abund)] <- 0
abund <- abund[,c(2:length(abund))]

#calculate FD using dbFD function, m set to 10 to avoid crach (int overload in QHULL)
FuncDiv <- dbFD(tr, abund, corr = c("lingoes"), m=10)

## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 49 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4209541

FDt <- as.data.frame(rbind(FuncDiv$FEve, FuncDiv$FDiv, FuncDiv$FDis, FuncDiv$RaoQ))
rownames(FDt) <- c("FEve", "FDiv", "FDis", "RaoQ")
colnames(FDt) <- c("A1", "A2", "A3", "A4", "A5", "A6", "A7")

FDt <- as.data.frame(t(round(FDt, digits=3)))

# create a nice table of FD indicies for the 7 management areas
formattable(FDt, list(    FEve = color_tile("transparent", "lightpink"),     FDiv = color_tile("transparent", "lightpink"),     FDis = color_tile("transparent", "lightpink"),    RaoQ = color_tile("transparent", "lightpink")))

	FEve	FDiv	FDis	RaoQ
A1	0.682	0.905	0.189	0.051
A2	0.637	0.849	0.263	0.079
A3	0.644	0.710	0.239	0.074
A4	0.614	0.745	0.200	0.048
A5	0.627	0.785	0.247	0.073
A6	0.665	0.650	0.179	0.048
A7	0.516	0.713	0.197	0.064

Rerunning the FD model excluding one trait at a time

Need to reduce ‘m’ to 9 to avoid crash:

Zero distance(s)Error in convhulln(tr.FRic, “FA”) : Received error code 2 from qhull. Qhull error: QH6114 qhull precision error: initial simplex is not convex. Distance=1.4e-16

Removing the MaxLenght trait from the model increased R2 by 0.047 compared to the model with all traits included. R2 was reduced when removing any of the other six traits from the model.

Conclusion

• Use the model without MaxLenght
• Removing Maxlength makes sence because the data included a few catches of sharks with very large MaxLengt (>2m) that skews the MaxLenght distribution. Excluding MaxLength creates an analysis around traits that are more common in distribution across species.
• Also, the sharks are rowing predators, not so site attached (but that also applies to some pelagic species, especially the large ones, so the same argument applies for these - excluding large body size reduces the importance of these non-site attached species)

# create data frame with original RaoQ and FRic
Raoq <- as.data.frame(c(FuncDiv$RaoQ, FuncDiv$qual.FRic))
colnames(Raoq) <- c("All")
rownames(Raoq) <- c("A1", "A2", "A3", "A4", "A5", "A6", "A7", "R2")


# select tr, but excluding one variable at a time
#tr[,c(1,!2,3,4,5,6,7)]
for (i in 1:7) { 
  FDout <- dbFD(tr[,c(setdiff(1:7, i))], abund, corr = c("lingoes"), m=9)
  raoq1 <- as.data.frame(c(FDout$RaoQ, FDout$qual.FRic))
  colnames(raoq1) <- colnames(tr)[i]
  Raoq <- cbind(Raoq, raoq1)
}

## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 37 PCoA axes (out of 46 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4684034 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 42 PCoA axes (out of 51 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3984771 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3663967 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3612455 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 48 PCoA axes (out of 57 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4029848 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.3538643 
## Species x species distance matrix was not Euclidean. Lingoes correction was applied. 
## FRic: Dimensionality reduction was required. The last 50 PCoA axes (out of 59 in total) were removed. 
## FRic: Quality of the reduced-space representation (based on corrected distance matrix) = 0.4132985

formattable(round(Raoq[8,], digits=3))

	All	MaxLength	Trophic.Level	Trophic.group	Water.column	Habitat	Gregariousness	diel
R2	0.421	0.468	0.398	0.366	0.361	0.403	0.354	0.413

formattable(round(Raoq[c(1:7),], digits=3), list(    MaxLength = color_tile("transparent", "lightpink"),     Trophic.Level = color_tile("transparent", "lightpink"),     Trophic.group = color_tile("transparent", "lightpink"),    Water.column = color_tile("transparent", "lightpink"), Habitat = color_tile("transparent", "lightpink"),    Gregariousness = color_tile("transparent", "lightpink"), diel = color_tile("transparent", "lightpink")))

	All	MaxLength	Trophic.Level	Trophic.group	Water.column	Habitat	Gregariousness	diel
A1	0.051	0.055	0.065	0.049	0.040	0.063	0.040	0.059
A2	0.079	0.102	0.103	0.059	0.059	0.093	0.070	0.089
A3	0.074	0.078	0.092	0.061	0.061	0.100	0.062	0.077
A4	0.048	0.063	0.063	0.050	0.031	0.046	0.039	0.057
A5	0.073	0.085	0.093	0.060	0.055	0.091	0.061	0.082
A6	0.048	0.050	0.058	0.036	0.042	0.064	0.046	0.051
A7	0.064	0.076	0.083	0.057	0.050	0.067	0.060	0.070

write.csv2(Raoq, "./figs/Raoq_table.csv")

---

Results table

This table combines CPUE,

st2 <- catch %>% group_by(ss) %>% summarise(first(id), first(Fhrs), first(gear), sum(weight))
colnames(st2) <- c("ss", "id", "Fhrs","gear", "weight")
st2$CPUEw <- st2$weight/st2$Fhrs

results.table <- data.frame(area=(1:7))
results.table <- cbind(results.table,(subset(st2, gear=="TB") %>% group_by(id) %>% summarise( mean(CPUEw)))[,2])
results.table <- cbind(results.table,(subset(st2, gear=="GN") %>% group_by(id) %>% summarise( mean(CPUEw)))[,2])
results.table <- cbind(results.table,(subset(catch.traits, gear=="TB") %>% group_by(id) %>% summarise( mean(TrophicLevel, na.rm = TRUE)))[,2])
results.table <- cbind(results.table,(subset(catch.traits, gear=="GN") %>% group_by(id) %>% summarise( mean(TrophicLevel, na.rm = TRUE)))[,2])
results.table<- cbind(results.table,(subset(catch.traits, gear=="TB") %>% group_by(id) %>% summarise( mean(as.integer(Gregariousness), na.rm = TRUE)))[,2])
results.table<- cbind(results.table,(subset(catch.traits, gear=="GN") %>% group_by(id) %>% summarise( mean(as.integer(Gregariousness), na.rm = TRUE)))[,2])
colnames(results.table) <- c("area", "TrapCPUE", "GillnetCPUE", "TrophicL_Traps", "TrophicL_Gillnets", "Greg_Traps", "Greg_Gillnets" )

TMean <- subset(catch.traits, gear=="TB" & TrophicLevel!="NA") %>% group_by(id, survey_m) %>% summarise(TM=mean(na.omit(TrophicLevel)))

results.table$SpeciesDensity <- sp.pr.area2$Species_Hrs
results.table$RaoQ <- Raoq$MaxLength[1:7]

formattable(round(results.table, digits = 3))

area	TrapCPUE	GillnetCPUE	TrophicL_Traps	TrophicL_Gillnets	Greg_Traps	Greg_Gillnets	SpeciesDensity	RaoQ
1	0.117	1.639	3.825	3.933	1.510	2.412	0.008	0.055
2	0.149	1.176	3.815	3.891	1.478	2.261	0.008	0.102
3	0.136	0.309	3.844	4.140	1.553	2.500	0.016	0.078
4	0.059	0.698	3.765	4.284	1.357	2.107	0.019	0.063
5	0.079	0.726	3.797	3.853	1.627	2.176	0.009	0.085
6	0.135	0.270	3.762	3.927	1.531	2.286	0.008	0.050
7	0.098	0.683	3.828	3.612	1.565	1.900	0.010	0.076

write.csv2(round(results.table, digits = 3), "./figs/results_table.csv")

---

Statistical analysis CPUE catches

Testing for normality of CPUE catches

Shapiro test for normality and Q-Q plots shows that data is non-normal and zero-inflated.

catch.df<-subset(catch1, number>0 & weight>0)

#' test for normality
shapiro.test(catch.df$CPUEw)

## 
##  Shapiro-Wilk normality test
## 
## data:  catch.df$CPUEw
## W = 0.5481, p-value < 2.2e-16

shapiro.test(catch.df$CPUEn)

## 
##  Shapiro-Wilk normality test
## 
## data:  catch.df$CPUEn
## W = 0.38853, p-value < 2.2e-16

#' Q-Q plots
library("qqplotr")

## 
## Attaching package: 'qqplotr'

## The following objects are masked from 'package:ggplot2':
## 
##     stat_qq_line, StatQqLine

gg <- ggplot(data = catch.df, mapping = aes(sample = CPUEw,  fill = factor(id))) + stat_qq_band() + stat_qq_line() + stat_qq_point(size=0.2) + facet_wrap(~ id, scales="free") +  labs(x = "CPUEw predicted", y = "CPUEw Observed") + ggtitle("Q-Q plots CPUEw")
gg

### Test for difference in CPUE between areas (non-parametric test) Both CPUEw and CPUEn are significantly different between the areas, but only CPUEn in the pairwise tests (nemenyi post-hoc) for areas 4:3 and 6:4. THis points to that Area 4 is the ‘odd one out’ in terms of CPUE along the Sudanses coast.

Significant difference for CPUEW between May 2013 survey and the Nov 2012 and Nov 2013 surveys, but not between the Nov12 and Nov13 surveys

#### NOT RUN BECAUSE ZERO-ADJUSTED GAM MODEL IS BETTER!

#' 1 - test for difference between areas
kruskal.test(catch.df$CPUEw, catch.df$id)
kruskal.test(catch.df$CPUEn, catch.df$id)

#' Post-Hoc tests of Kruskal - Wallis results
library("PMCMR", lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")
posthoc.kruskal.nemenyi.test(catch.df$CPUEw, catch.df$id, "Chisq")
posthoc.kruskal.nemenyi.test(catch.df$CPUEn, catch.df$id, "Chisq")
#' significant difference in CPUEn between area 6 and area 4

#' 2 test for differences between surveys
kruskal.test(catch.df$CPUEw, catch.df$survey)
kruskal.test(catch.df$CPUEn, catch.df$survey)

posthoc.kruskal.nemenyi.test(catch.df$CPUEw, catch.df$survey, "Chisq")
posthoc.kruskal.nemenyi.test(catch.df$CPUEn, catch.df$survey, "Chisq")

Zero-inflated model (better approach)

This method is more appropriate than the Kruskal-Wallis approach above as it evaluates the effects of all factors at the same time, instead of evaluating the effects of area separate from survey.

Developed a model for CPUEw including depth, area, survey and gear, which performed better than a simpler model excluding gear.

Depth is not significant

Significant factors:

• Area (1, 5, 6 and 7)
• Surveys (all)
• Gear (all gear types)

#' Fitting GAM Zero Inflated model to catch data

# select only variables used in model
c2 <- catch[,c(1,14,22,25,28)]
c2 <- subset(c2, gear=="TB" | gear=="GN" |gear=="HL")

# CPUEw - GAMLSS ZAGA model incluing survey, area, depth & gear 
mod<-gamlss(CPUEw~depth+factor(id)+factor(survey)+factor(gear),family=ZAGA, data=c2) 

## GAMLSS-RS iteration 1: Global Deviance = 279.0768 
## GAMLSS-RS iteration 2: Global Deviance = 279.0767

summary(mod)

## ******************************************************************
## Family:  c("ZAGA", "Zero Adjusted GA") 
## 
## Call:  gamlss(formula = CPUEw ~ depth + factor(id) + factor(survey) +  
##     factor(gear), family = ZAGA, data = c2) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           -1.183870   0.181345  -6.528 9.87e-11 ***
## depth                  0.003163   0.002327   1.359 0.174344    
## factor(id)2           -0.066567   0.149632  -0.445 0.656495    
## factor(id)3           -0.076983   0.202036  -0.381 0.703246    
## factor(id)4           -0.068939   0.202976  -0.340 0.734189    
## factor(id)5           -0.513933   0.167731  -3.064 0.002233 ** 
## factor(id)6           -0.401688   0.162065  -2.479 0.013331 *  
## factor(id)7           -0.585171   0.164508  -3.557 0.000390 ***
## factor(survey)2013002 -0.162780   0.095859  -1.698 0.089748 .  
## factor(survey)2013005 -0.224448   0.111011  -2.022 0.043418 *  
## factor(gear)HL         0.775526   0.194477   3.988 7.08e-05 ***
## factor(gear)TB        -0.415457   0.108136  -3.842 0.000129 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.08571    0.02157   3.974 7.49e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.74653    0.06204  -12.03   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1191 
## Degrees of Freedom for the fit:  14
##       Residual Deg. of Freedom:  1177 
##                       at cycle:  2 
##  
## Global Deviance:     279.0767 
##             AIC:     307.0767 
##             SBC:     378.2324 
## ******************************************************************